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PARTICLE-IN-CELL SIMULATION OF ELECTRON ACCELERATION IN SOLAR CORONAL JETS 
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ABSTRACT 

We investigate electron acceleration resulting from 3D magnetic reconnection between an emerging, 
twisted magnetic flux rope and a pre-existing weak, open magnetic field. We first follow the rise 
of an unstable, twisted flux tube with a resistive MHD simulation where the numerical resolution is 
enhanced by using fixed mesh refinement. As in previous MHD investigations of similar situations, the 
rise of the flux tube into the pre-existing inclined coronal magnetic field results in the formation of a 
solar coronal jet. A snapshot of the MHD model is then used as an initial and boundary condition for a 
particle-in-cell simulation, using up to half a billion cells and over 20 billion charged particles. Particle 
acceleration occurs mainly in the reconnection current sheet, with accelerated electrons displaying a 
power law in the energy probability distribution with an index of around -1.5. The main acceleration 
mechanism is a systematic electric field, striving to maintaining the electric current in the current 
sheet against losses caused by electrons not being able to stay in the current sheet for more than a 
few seconds at a time. 

Subject headings: acceleration of particles — Sun: corona — Sun: magnetic topology 



1. INTRODUCTION 

Solar jets have been shown to be triggered by mag- 
netic reconnection, similarly to solar flares, while their 
released energy is much below what is set free in a 
medium-sized flare event, and the timescales are usu- 
ally shorter. Nevertheless their high frequency of oc- 
currence makes them a significant contributor to the so- 
lar ejecta, particularly the solar wind originating from 
coronal holes. Solar jets feature upflow velocities of 
more than 150 km s" 1 (jSavcheva et al.ll2007llChifor et al.l 
2008) and are observed at EUV d own to X-ray wave - 
lengths primarily in corona l holes (|Kamio eTaI][2003), 
but also in ac tive regions (|Chifor et al.l |2008T h Using 
RHESSI data, iKrucker et al.l (|2007l 120111 ) have further 
investigated the electron impact regions of solar jets as 
a result of interchange reconnection. 

There have been several studies in the past employ- 
ing fully 3D kinetic models to study particle accelera- 
tion. Two approaches are most popular: particle-in-cell 
(PIC) simulations and test particle simulations. Their 
main difference is the back-reaction from the particles 
onto the fields, which is only taken into account in the 
former method, while it is assumed to be negligible in 
the latter one, justified by using a low number of test 
particles for such simulations. These kind of simulations 
provide much information on particle tr ajectories and fa- 
vored locations of particle acceleration (iTurkmani et al.l 
[2^051 [20CjllDalla fe Browningll2005[ I2006L 120081) . Never- 
theless there are severe limitations to test particle sim- 
ulations. These together with their consequenc e s hav e 
in detail been discussed in Ros dahl fc Galsgaardl (|2010( ). 
On the other hand, the approach using self-consistently 
evolving fields such as in PIC codes is subject to several 
resolution constraints, reducing the possible physical box 
size that can be simulated to far below length scales of 
solar jets. In order to bypass these limitations, we made 
use of modifications of the constants of nature. Thereby 
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we are able to present results from fully 3D PIC simu- 
lations of a solar jet, using essen t ially the same initial 
setup as in iRosdahl fe G alsgaardl (|2010h . but with self- 
consistently evolving fields. Such mo difications have pre - 
vio usly to some extent been use d bv lDrake et al.l (2006) 
and lSiverskv fe Zharkov a (2009). A comparison of simu- 
lations using different modifications ex ceeds the scope 
of thi s Letter. We therefore refer to I Baumann et al.l 
(2012) for a description and analysis of the physical con- 
sequences of modifications. For the present study we use 
a choice of modifications as close as we can get to reality 
with the currently available computational resources. 

Section [2] provides an overview of the experiment and 
describes the MHD and PIC simulations as well as their 
intercoupling. In Section [3] the results are presented and 
discussed. Finally, in Section |4] conclusions are drawn. 

2. SIMULATIONS 

The solar jet experiment at hand starts out with a 
fully 3D resistive compressible MHD simulation of a 
twisted emerging flux rope, initially positioned 1.7 Mm 
below the photosphe r e, sim ilar to the setup used by 
iMoreno-Insertis et al. (2008). A constant magnetic field 
of 3.3 G is imposed on the entire computational box, in- 
clined 65° in the yz-plane. The maximum magnetic field 
strength of the flux rope is slightly higher than 1000 G 
and hence much larger than the background magnetic 
field. The atmosphere is initially in hydrostatic equilib- 
rium wi th a ID atmosp h eric p rofile similar to the one 
used in lArchontis et al.l (|2005f ): the sub-photospheric 
temperature at the bottom is 5.5xl0 4 K, with a maxi- 
mum mass density p of about 9xl0 _6 gcm -3 at a depth 
of 3.7 Mm below the surface. The "chromosphere" has a 
constant temperature of around 5600 K and the corona 
starts out with T = 2.2 x 10 6 K and p = 6xl0" 16 gcm" 3 , 
as illustrated in Figure[TJ 

The simul ations are performed us i ng th e Stagger MHD 
code, as in IMoreno-Insertis et al.l ([2008), assuming an 
ideal gas law and neither taking heat conduction nor 
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Fig. 1. — ID atmospheric profile for the MHD simulations. The 
dashed vertical line shows the lower cut in height for the sub- 
domain used for the PIC simulations. 

radiative cooling into account. Viscosity and resistiv- 
ity are locally defined, depending mainly on the velocity 
divergence, in order to provide a suitable (to maintain 
code stability), but minimal amount of dissipation. A 
mesh with 512 3 cells covers a box with physical extents 
of 33.8 x 38.1 x 32.5Mm. The mesh is non-uniform in all 
directions, with a minimal mesh spacing of (x m i n , y m im 
Zmm) = (0.034, 0.034, 0.030) Mm around the reconnection 
region, and with a mesh spacing less than 10% larger 
than that in a region of size 12.5 x 10.7 x 5.1 Mm. The 
z-coordinate is the direction normal to the solar surface. 

Despite the slightly denser coron a compared to the pre- 
viousl y performed simulations by iMoreno-Insertis et al.l 
(2008), the overall evolution of the experiment is the 
same: the twisted flux tube is made buoyantly unstable 
by applying a density perturbation at its center, which 
causes it to rise up into the much rarer corona, against 
the Lorentz force fr om the bending t ube, w hile expand- 
ing as described in lArchontis et al.l (|2004D . Above the 
photosphere, the expansions of the magnetic flux rope, 
now due to the high magnetic pressure, continue rapidly 
as more and more flux reaches coronal heights. At the 
same time the corona is locally pushed upward where 
pla sma as well as magn etic flux emerges from below 
(cf. lArchontis et ail 120051 ). The interaction between the 
two magnetic field domains defined by the initial coronal 
magnetic field and the flux rope leads to a destabiliza- 
tion of the field configuration and causes the formation 
of a thin dome-shaped current sheet where the magnetic 
field lines are most inc l ined r elative to each other. As in 
IMoreno-Insertis et al.l ([2008), most magnetic field lines 
of the two domains end up being almost anti-parallel at 
their first encounter, which makes their interaction max- 
imally powerful. The sheet is subject to ohmic dissipa- 
tion, causing it to reach temperatures as high as about 
8 x 10 6 K. Reconnection gradually occurs between various 
field lines along the twisted tube and the ambient coronal 
magnetic field. Due to this restructuring of the magnetic 
field, sev eral distinct flux doma i ns for m in the corona, as 
shown in lMoreno-Insertis et al.l ((20081 see their Figure 2). 
A hot plasma jet pair emerges from the high tempera- 
ture and gas pressure region of the reconnection region, 
propagating sideward together with the expelled recon- 
nected magnetic field lines. The plasma in the jets is fed 
to the reconnection site from both sides of the current 
sheet; the region below supplies dense and cold photo- 



spheric plasma, while the plasma coming from above is 
much hotter and rarer coronal plasma. 

After a large fraction of the dense plasma enclosed in 
the flux rope has been reconnected and ejected in the 
form of plasma jets, as well as drained off along the mag- 
netic field lines due to gravity acting on the heavy sub- 
photospheric gas, we initialize the PIC experiment with 
a cut-out of size 22 x 22 x 22 Mm from the MHD data 
set. The reconnection process continues in the PIC sim- 
ulation expelling coronal low-density plasma and recon- 
necting field lines from both connectivity domains. 

The PIC sim ulations are performed using the Photon- 
Plasma code (|Haugb0lld [2001 iHededall I2005T ). which 
solves the Maxwell equations together with the relativis- 
ts equation of motion for charged particles. We fix the 
magnetic fields at the boundaries to the values given by 
the MHD data set, and lea ve the boundaries ope n for 
particles to exit and enter (Haugb0llc et al. 2012). To 
initialize the electric field in the PIC simulation, only 
the advective electric field (— u x B, where u is the bulk 
speed) is passed on. The particles are initially given a 
random thermal velocity drawn from a Maxwellian distri- 
bution, plus the bulk velocity from the MHD simulation. 
Electron velocities consist of additionally the velocity due 
to the initial electric current: 



vj = 



-(VxB), 



(1) 



where the magnetic field B and the density n are pro- 
vided by the MHD snapshot data set. The jet velocity 
of this specific MHD data set is at that point in time on 
the order of 400-800 km s _1 . These jets are dominated 
by the thermal motion in this high-temperature and low 
magnetic field region. 

We conducted several PIC runs with grid dimensions 
of 400 3 and 800 3 cells on a uniform grid with cell sizes 
of 55 km and 27.5 km respectively, in each case with 20 
particles per species (protons and electrons) per cell, sim- 
ulating up to 7.5 solar seconds. 

To minimize computational constraints, the MHD 
snapshot is cut at 1.1 Mm below the bottom of the 
corona, hence in the transition region, shown as a vertical 
dashed line in Figure[T] This limits the density span to 
a factor of 4xl0 4 , small compared to the span of about 
2xl0 10 covered by the MHD simulation, but still large 
for a PIC simulation. Additionally, the temperature is 
reduced by a factor of four in order to avoid drowning 
the high-energy tail in the Maxwellian distribution. 

To make the plasma micro-scales marginally resolv- 
able, the charge per particle is reduced, and to ease 
the time step constraint from the propagation of electro- 
m agnetic waves the speed of light is reduced, as explained 
in iBaumann et al.l ([20121 ) . The electron skin depth in 
the current sheet is resolved with about 5-10 grid cells. 
The electron gyroradius varies between a fraction of a 
cell in the flux tube interior to many cells inside the 
current sheet — because of the near cancellation of oppo- 
sitely directed magnetic fields there are several dynami- 
cally evolving null points inside the current sheet. 

In addition to these stratified atmosphere simula- 
tion runs, control runs with a constant density of 
3xl0 -15 gem -3 were also performed. Both types of runs 
show similar overall results. 
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Fig. 2. — Magnetic field lines (red/black) with (a) the diffusive electric field E ros = rrj in the PIC cut-out of the MHD snapshot data set 
(purple/gray volume) and the charge density plane at the bottom of the box — note the low density inside the flux rope, as explained in the 
text — and (b) the Eu field 0.5s after start of the PIC simulation in run 400 3 (purple/gray), together with the current density plane at the 
bottom of the box. 

3. RESULTS AND DISCUSSION 

During the MHD flux emergence simulation a diffusive 
electric field E ros = 77 j , where rj is the resistivity and j 
the electric current density, builds up in the reconnection 
region, approximately cospatial with the current sheet. 
The parallel (in relation to the magnetic field) electric 
field £j| is part of the diffusive (non-ideal) electric field, 
and provides information on the rate of reconnection as 
well as on favored regions for particle acceleration. The 
diffusive component of the electric field is, unlike the ad- 
vective electric field (— u x B), not inherited by the PIC 
code, but builds up self-consistently. Figure[3] compares 
the location of the diffusive electric field component for 
the chosen snapshot of the MHD simulations with the Eu 
field of the PIC simulation 0.5 s after start. The parallel 
electric fields reach in general much higher magnitudes in 
the MHD simulations compared to the PIC simulations. 

Eu is the most efficient particle accelerator, since its 
force acts on the particles without being affected by the 
perpendicular particle gyromotion. Its maximum in the 
PIC simulations is located inside the current sheet, equiv- 
alent to the diffusive electric field E rcs in the MHD sim- 
ulation. Accelerated electrons (see Figure[3]) are located 
in the plasma outflow parts of the reconnection region. 
The electron bulk velocity in the jet is on the order of 
2000 km s -1 . The proton bulk jet flow on the other hand 
is only about 270 km s -1 , which difference to the elec- 
tron bulk speed defines the electric current required by 
the magnetic field configuration. The lower plane of the 
PIC simulation visualization in Figure^] (right) shows the 
electric current density. At the bottom center of this fig- 
ure resides the flux rope, whose twisted field lines are 
indicated by the strongest electric current pattern. Addi- 
tionally, to the left of the flux rope signature, the current 
sheet features a turbulent structure. This is the result of 
the plasma transport caused by reconnection. In addi- 




FlG. 3. — Electrons that win energy over a time interval of 2.5 s, 
about 4.5 s after the start of the 400 simulation run together with 
the magnetic field lines (white) and the electric current density 
as contour plot (blue-green-red/black-gray-white for increasing 
current density) in a j/2-plane. Particles with velocities directed 
upward are colored purple /light gray, while downward moving elec- 
trons are colored green/gray. 

tion, we observe fast up and down flowing plasma, as can 
be seen in Figure[3l in which upward moving electrons are 
shown in purple/light gray and downward moving elec- 
trons are shown in green/gray. 

By tracing particles that win energy over a time pe- 
riod of a second, it becomes clear, see FigurelU that the 
acceleration mechanism is a systematic DC electric field, 
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Fig. 4. — Four random electrons traced in run 800^ during Is. Their projected positions are plotted in the slices to the right, together 
with the electric current density (raised to the power 0.5 to enhance fine structures) in yz- and jyx-planes. The black dashed lines in the 
right images show the cut in the respective direction for the other image. Additionally, the gyroradius r g of the particles, their cosine of 
the pitch angle cos(B,v), and the electric field in the direction of motion, E\\v, are plotted. 



particularly present in the proximity of the current sheet, 
and mainly directed parallel to the magnetic field. How- 
ever, since the magnetic field is very weak in the current 
sheet area particles are almost decoupled from the mag- 
netic field lines, leading to pitch angle fluctuations seen 
in the third panel to the left in FigureHJ 

We only write out every 25th particle data, hence in 
Figures|3] and [4] one particle represents in reality 25 of 
its type. Note that in Figure^] the electric field fluctua- 
tions in the second panel on the left are to a large extent 
Monte Carlo noise, due to low numbers of particles per 
cell. This has been verified with a control experiment 
using twice the number of particles per cell. Further, 
the background electric current density image is not ex- 
actly at the particle positions, which results in a slight 
projection offset. 

Only a limited number of particles happen to be close 
enough to the current sheet to be efficiently accelerated 
by E\\ , thus being able to contribute to the electric cur- 
rent required by the magnetic field topology. After ex- 
periencing a certain amount of acceleration, they get ex- 
pelled from the current sheet, as they follow magnetic 
field lines which leave the current sheet region. The par- 
ticles lost from the current sheet are replaced by new par- 
ticles, which again need to be accelerated. Since the mag- 
netic field in the reconnection region is very low, due to 
the chosen orientation of the initial background magnetic 



field geometry, it is easy for electrons to get misguided, 
as they are no longer tightly attached to their magnetic 
field lines and therefore may encounter different electric 
field structures on their large gyroradius trajectories. 

The systematic parallel electric field building up in 
the current sheet area is capable of accelerating par- 
ticles to non-thermal velocities. Figure[5] presents the 
energy histogram of electrons located in a cut-out of 
19.0 x 13.2 x 6.5 Mm around the reconnection region. 
The initial energy distribution is shown by the dashed 
line. It is primarily a superposition of the two different 
plasma inflow domains of the reconnection region passed 
to the PIC code through the MHD bulk velocity: one of 
coronal origin, hence at higher temperatures and lower 
densities, and another one from plasma emerging from 
the flux rope, hence at lower temperatures and higher 
densities. To this the drift speed from particles in the 
current sheet is added, as defined by Equation (fTJ) . The 
distribution of probability for a particle to stay in the 
current sheet defines a power law in the high-energy tail 
of the Maxwell-Boltzmann distribution, typically featur- 
ing an index of around -0.5 on a logarithmic scale. This 
is the case for electrons and protons, as the acceleration 
mechanism inside the current sheet is for both species the 
same. A series of different simulations indicates, that 
the influence of the physical resolution on the power- 
law index as well as on the upper energy limit is mi- 
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Fig. 5. — Left: normalized electron probability distribution at t = 4s, for particles in a cut-out of size 19.0 X 13.2 X 6.5 Mm around the 
reconnection region for a 400 3 and a comparable 800 3 simulation and a zoom-in. The black dashed curve illustrates the initial particle 
distribution. The vertical dotted lines mark the area for which power laws are fitted, whose indices are annotated. The lower energy limit 
coincides with the lowest considered electron energy for the image to the right. Right: the sum of the deposited electron energy in the 
lowest Mm of the computational box collected over a time period of 7s (red, low; blue, high). The viewing angle is as in Figuref2] but seen 
from the top. 

ulations. 

A strong correlation is found between a slowly evolving 
DC electric field located inside the current sheet and the 
location of the accelerated particles (electrons and pro- 
tons). The magnetic field is weak and chaotic inside the 
current sheet, from where most of the particles that are 
accelerated are quickly lost, only to be replaced by new 
particles, which again need to be accelerated. This repre- 
sents the dissipation mechanism at work. The systematic 
electric field required to constantly accelerate new par- 
ticles is, in effect, a dissipative ("resistive", non-ideal) 
electric field, sustained even though the plasma particles 
in this experiment are collisionless. This physical process 
leads to an energy probability distribution of accelerated 
particles, featuring a power-law tail at the high-energy 
side of the initial Maxwell-Boltzmann distribution, with 
a power-law index of about —1.5. The power-law index 
and the invariant energy cut-off show numerical conver- 
gence of the simulations. 

The impact region of energetic electrons on the lower 
boundary of the box is comparable to observations of 
si milar events and mag netic field geometries, described 
bv lKrucker eTaTl (|201lf ). 

In the future, with increasing available computational 
resources, we hope to be able to resolve turbulent re- 
gions sufficiently, to enable a study of stochastic particle 
acceleration, expected to be able to accelerate particles 
to much higher energies. 
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nuscule (see Figure[5|), since the current per unit width 
of the current sheet is set by the magnetic field geome- 
try. This independence also demonstrates numerical con- 
vergence. The tail slopes converge very rapidly toward 
these power-law indices, since the electron acceleration 
happens impulsively. Part of the accelerated electrons 
escape along open magnetic field lines (see Figure[3|) into 
interplanetary space and hence contribute to the solar 
wind. Downward moving energetic electrons impact the 
much denser chromosphere, which could on the Sun lead 
to bremsstrahlung emission. In the case of high-energy 
electrons, this may cause hard X-ray signatures of the 
type observed by RHESSI. Figure[5] shows that most en- 
ergy is deposited at the magnetic field line footpoints of 
the coronal loops. Obser vational hard X - ray si gnatures 
of solar jets described in iKrucker et al.l (|201 If ) show a 
three-footpoint pattern (see, e.g., event 7 in Figure 3). 
Due to the choice of our flux rope being 90° inclined 
relative to the uniform background magnetic field, we 
rather get a 3D structure of footpoints, represented by 
three lines of electron impact regions in Figure[SJ Ad- 
ditionally, because of the magnetic field line orientation 
and twist, an impact region at the intersection of the 
flux tube and the photosphere develops. (The impact 
line to the right is slightly cut because of the choice of 
the computational box size.) 

4. CONCLUSIONS 

On the basis of an MHD jet experimen t , simi lar to the 
one conducted by iMoreno-Insertis et al.l (|2008|) . but us- 
ing stretched meshes to obtain higher spatial resolution, 
we have used PIC simulations to study the acceleration 
of charged particles in the 3D reconnection region of a 
solar coronal jet. This is the first fully 3D kinetic model 
employing self-consistent fields to investigate particle ac- 
celeration in the context of solar jets. It uses a new 
concept of combining macroscopic with microscopic sim- 
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